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ABSTRACT 

We study the interaction between tides and convection in astrophysical bodies by 
analysing the effect of a homogeneous oscillatory shear on a fluid flow. This model can 
be taken to represent the interaction between a large-scale periodic tidal deformation 
and a smaller-scale convective motion. We first consider analytically the limit in which 
the shear is of low amplitude and the oscillation period is short compared to the 
timescales of the unperturbed flow. In this limit there is a viscoelastic response and 
we obtain expressions for the effective elastic modulus and viscosity coefficient. The 
effective viscosity is inversely proportional to the square of the oscillation frequency, 
with a coefficient that can be positive, negative or zero depending on the properties 
of the unperturbed flow. We also carry out direct numerical simulations of Boussinesq 
convection in an oscillatory shearing box and measure the time-dependent Reynolds 
stress. The results indicate that the effective viscosity of turbulent convection falls 
rapidly as the oscillation frequency is increased, attaining small negative values in 
the cases we have examined, although significant uncertainties remain because of the 
turbulent noise. We discuss the implications of this analysis for astrophysical tides. 
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1 INTRODUCTION 

Tidal interactions determine the orbital and spin evolution 
of astrophysical bodies when they orbit sufficiently close to 
one another. Applications include close binary stars, extra- 
solar planetary systems and the satellites of solar-system 
planets. 

In many cases of interest the tidally forced bodies are 
fully or partly convective. In order to determine the rate of 
tidal evolution it is therefore necessary to study the interac- 
tion between tides and convection. The response of a fluid 
body to tidal forcing generally consists of two components: 
one is non-wavelike and of large scale, while the other in- 
volves internal waves of smaller scale. At the simplest level 
the convection could be thought of as providing an effective 
viscosity that damps the non-wavelike tidal disturbance and 
provide s a ph ase shift in the response. As in the theory of 
iDarwinl (18801 ). the tidal torque and the rate of tidal evolu- 
tion are then directly proportional to this effective viscosity. 
Convection may also play an important role in dissipating 
inertial waves, which constitute the low-frequency wavelike 
respons e of rotating convecti ve zones of stars and giant plan- 
ets (e.g. lOgilvie fc Linll2004h . 

As pointed out by IZahnl (1 196ot ) and 



iGoldreich fc Nicholson! |l977]), the effective viscosity 
estimated from mixing-length theory ought to be reduced 
when the period of the tidal disturbance is short compared 
to the characteristic timescale of the convective motion. 
The suppression factor has been the subject of much 
debate, informed by observational constraints such as the 
appar ent rate of circularization o f the orbits of close binary 
stars ( Meibom fe Mathieul 120051 . and references therein). 
IZahnl (fl966h suggested a less severe suppression, such that 
the effective viscosity is proportional to the oscillation 
period for short periods, which is in better ag re ement 

IGoldreich fc Nicholson! (|l977T ) 



with these observations. 



and IGoldreich fc Keelevl (| 19771 ) considered the multiscale 
nature of the turbulent flow, and argued that the dominant 
contribution to the effective viscosity at short periods 
comes from eddies with a convective timescale comparable 
to the oscillation period. For a Kolmogorov spectrum, this 
argument gives a more powerful suppression, such that 
the effective viscosity is proportional to the square of the 
oscillation period for short periods. 

While these authors relied on si mple physical argument s 
and order-of- magnitude estimates. iGoodman fc Ohl (|l997l h 
who also provide a clear review of the controversy, intro- 
duced a more formal procedure for determining the effec- 
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tive viscosity of a convective flow, by considering the effect 
of a homogeneous oscillatory strain on a turbulent velocity 
field. Analytical progress was limited by the complexity of 
the equations. While their method relies on an expansion 
in powers of the ratio of the oscillation period to the con- 
vective timescale, it leads to a result in which the dominant 
contribution comes from eddies for which these timescales 
are similar. Their ar gumen t does , however, appear to rule 
out the hypothesis of IZahnl (|l966T l. 

More recently, numerical simulations o f convection 
have been brought to b e ar on this question. iPenev et al.l 
(120071) and iPeney et al.l (|2009l ) applied the procedure of 
iGoodman &i Oh to the velocity fields obtained in 

numerical simulations of convection in a deep layer, to de- 
duce the effective viscosity as a funct i on of the tidal fre- 
quency. In lPenev. Barranco fc Sasselovl l|2009h an oscillatory 
forcing was introduced directly into the convection simula- 
tion and the effective viscosity was estimated by measuring 
the work done by this force. The results of these s tudies 
sugges t that something closer to the prescription of IZahnl 
(1966) may be appropriate, not for the reasons originally 
suggested, but possibly because the power spectrum of the 
convectio n is less steep than the Kolmo gorov spectrum as- 
sumed bv lGoldreich fc Nicholson! l|l977h . 

While the results of Penev and collaborators are of 
considerable interest, the uncertainties in these calculations 
have not been quantified. Owing to the importance of this 
problem, we are motivated to examine it from an indepen- 
dent viewpoint. We have devised a theoretical framework 
that allows us to study statistically homogeneous flows; 
while this local viewpoint omits important global aspects 
of stellar or planetary convection, it allows us to resolve 
convection at high Rayleigh numbers and to quantify the 
uncertainties due to turbulent noise. We also present com- 
plementary analytical results which shed new light on pre- 
vious theoretical discussions. 

The plan of this paper is as follows. In Section [2] we 
introduce the homogeneous sheared coordinate system used 
in our analytical and numerical calculations, and discuss the 
equations of fluid dynamics in this system. In Section [3] wc 
derive the response of a fluid flow to oscillatory shear in the 
limits of low amplitude and high frequency, and calculate 
the effective elasticity and viscosity in various cases. Direct 
numerical simulations of convection in an oscillatory shear- 
ing box are reported and interpreted in Section [4] followed 
by a summary and discussion of the results. 



2 THE OSCILLATORY SHEARING BOX 
2.1 Motivation 

At the present time it is not practical to compute a global 
model of an astrophysical body with turbulent convection 
and to measure its tidal response through the direct appli- 
cation of a time-dependent gravitational potential. As noted 
in the introduction, the response of a fluid body to tidal forc- 
ing generally consists of a large-scale non-wavelike motion 
together with some internal waves. The non-wavelike tide 
can be computed without difficulty in a linear approxima- 
tion and consists of a time-dependent spheroidal deforma- 
tion of the body. From the local perspective of the smaller- 
scale convective motion, this tidal flow appears as a spatially 



homogeneous, approximately incompressible motion that is 
oscillatory in time. We are interested, therefore, in deter- 
mining the response of the convective motion to such a de- 
formation and the additional stresses that result from this 
interaction. 

A general, three-dimensional, spatially homogeneous, 
incompressible deformation can be represented by a traceless 
velocity gradient tensor Vu that is independent of position. 
In this paper we consider a time-dependent simple shear 
of the form u y oc x. Using a linear combination of simple 
shears with different orientations it is possible to compose 
an arbitrary velocity gradient tensor that is traceless (see 
AppendixfX]). Provided that we are working a linear regime, 
therefore, it is sufficient to measure the response to a time- 
dependent simple shear. 

2.2 Sheared coordinates 

We initially consider an incompressible fluid in two dimen- 
sions, satisfying the Navier-Stokes equations 

(d t + u x d x + u y d y )u x = -d x p + v(d 2 x + d 2 )u x + f x , (1) 

(8 t + U X d X + Uydy)Uy = ^OyP + u(6^ + 8y)Uy + fy, (2) 



B X U X + dyUy — 0, 



(3) 



where u is the velocity, p is the pressure divided by the 
density, v is the kinematic viscosity and / is a body force 
per unit mass. 

We suppose that the system is subject to a homoge- 
neous time-dependent shear. We define the sheared coordi- 
nates 



x — x, y = y — a(t)x, t = t. 



(4) 



The dimensionless quantity a is the shear, and d = da/dt 
is the shear rate. The Jacobian determinant of the trans- 
formation is unity. Partial derivatives transform according 
to 

dx = d' x - ad'y, d y = d' y , d t = d' t - axd' y . (5) 
We write 

u x —v x , u y —Vy + ax, (6) 

so that v represents the velocity relative to the shearing mo- 
tion. We continue to refer vector components to the original 
Cartesian basis. 

With these substitutions, the Navier-Stokes equations 
become 

[d' t + v x (d x -ad y ) + v y d' y ]v x = -(d' x -ad' y )p+vAv x + f x , (7) 
[d' t + v x (d' x — ad'y) + Vydy]v y + ax + av x = — d' y p 

+ VAVy + fy, 

(d' x - ad' y )v x + d' y Vy = 0, (8) 
involving the Laplacian operator in unsheared coordinates, 

a = dl + dl 



= {d' x - ad'y) 2 + d. 



i \2 , q/2 

y 

A' — 2ad' x d' y + a 2 d'y. 
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These equations are spatially homogeneous (not involv- 
ing x' or y explicitly) except for the term ax' , which arises 
because a spatially inhomogeneous force is required to main- 
tain an accelerating shear. By writing 

f*=fx, fy=ax' + f y , (9) 

i.e. by subtracting the inhomogeneous force, we restore spa- 
tial homogeneity to the equations. In this representation, / 
is the body force that, if necessary, drives the (possibly tur- 
bulent) motion whose response we wish to measure, while 
ax' e y is the force that drives the imposed large-scale shear. 

The property of spatial homogeneity results from the 
translational invariance of the homogeneously sheared sys- 
tem, given that a Galilean transformation easily removes 
the velocity shift associated with a translation in the x di- 
rection. It means that statistically homogeneous turbulence, 
for example, is possible in this system. It also means that, for 
computational or analytical purposes, we can impose peri- 
odic boundary conditions on v and p in sheared coordinates. 
When employing a spectral method, these variables can be 
expanded in Fourier series in sheared coordinates, using ba- 
sis functions exp(ifc^a;' + ik' y y'). 

An alternative approach, which leads to identical re- 
sults, is to solve the original equations in unsheared coor- 
dinates but to apply modified periodic boundary conditions 
of the form 

p(L x , y, t) = p(0, (y - a(t)L x ) mod L y ,t), (10) 

p(x,L y ,t) =p(x,0,t), (11) 

and similarly for v (not u), where L x and L y are the di- 
mensions of the shearing box. The appropriate basis is then 
one composed of shearing waves exp(ik x (t)x + ik y y), having 
(quantized) time-dependent wavevectors with components 
k x (t) = k' x — a(t)k' y and k y = k' y . With these definitions, 
exp(ik x (t)x + ik y y) = exp(ik' x x' + ik' y y'). 

Note that the oscillatory shear in our model is not im- 
posed by the boundary conditions as such, but is driven by 
the inhomogeneous body force described above. 

Our model therefore provides a self-consistent local repre- 
sentation of any turbulent flow subject to an oscillatory de- 
formation due to an external body force. 

Similar considerations apply to the Boussinesq equa- 
tions in three dimensions, or indeed to the equations of 
compressible convection. In the Boussinesq case the basic 
equations in unsheared coordinates are 

(d t + Ujdj)u-i = —dtp + vAm + bd + f i% (12) 
(d t + u l d l )b + N 2 u z e z = KAb, (13) 

d z u, = 0, (14) 

where 6 is a buoyancy variable (proportional to the Eulerian 
entropy perturbation multiplied by the gravitational acceler- 
ation), e is a unit vector in the direction opposite to gravity, 
TV 2 is the square of the buoyancy frequency (negative in the 
convectively unstable case) and k is the thermal diffusiv- 
ity. (In the case of convection no additional body force / is 
required to drive the motion.) 

The sta n dard s hearing box, as e m ployed by, 
e.g., iRogallol (|l98lft . IWisdom fc Tremaind \l98$ ) and 



lHawlev. Gammie fc Balbusl (|l995h . corresponds to setting 
a(t) oc t and, if appropriate, adding rotation to the box. 
In this case the shear is inexorable and the shearing 
coordinates must be remapped periodically for the purposes 
of numerical simulation. 

In our case we consider a(t) to oscillate sinusoidally and 
no remapping is required. As the frequency of the oscillation 
is varied, we can choose either to keep the amplitude of a 
the same, which means that the maximum angle of the shear 
is fixed, or to scale the amplitude of a inversely with the 
frequency, which means that the maximum angular velocity 
of the shear is fixed. For the most part we are interested in 
the regime in which the shear is very small (in either sense), 
but if it is too small then its effects cannot be measured 
reliably in a direct numerical simulation. 

The quantity to be measured that is of greatest interest 
is the shear stress associated with the fluid motions, i.e. the 
Reynolds stress component 

- R xy = —{v x v v }, (15) 

where the angle brackets denote a suitable averaging opera- 
tion. (The physical Reynolds stress has an additional factor 
of the density of the fluid.) It is this stress that exchanges 
energy with the large-scale shear. The energy equation in 
sheared coordinates in the absence of buoyancy forces is 

d' t {\viVi) + diFi = -av x v y - v\ V x v\ 2 + fiVi, (16) 

where F is an appropriate energy flux density. If the 
Reynolds stress —R xy is positively correlated with the shear 
rate a then, like a viscous stress, it extracts energy from the 
large-scale shear and imparts it to smaller-scale motions. 
However, there is no reason in principle why the correlation 
should not be negative, in which case work is done on the 
large-scale shear at the expense of either the body force or 
the decaying energy of the smaller-scale flow. 



3 HIGH-FREQUENCY RESPONSE 
3.1 Linearized equations 

In this section we present analytical results on the response 
of fluid motions to high-frequency shear of low amplitude. 
We omit buoyancy forces, which are probably inessential for 
this discussion, and consider the Navier-Stokes equations in 
sheared coordinates, 

[d' t + v ] (d' J - aS ]1 d' y )]v i + av x 8 i2 = ~{d[ - aSnd^p 
+ v(d'j - a8jid' y )(d'j - adjid' y )vi + 

(3- - a5 ll d' y )v i = 0, (17) 

where Sij is the Kronecker delta and the subscripts 1 and 2 
refer to the x and y directions involved in the shear. We now 
carry out a linearization in the shear amplitude. We consider 
a basic flow that exists in the absence of shear, satisfying the 
equations 

(d' t + Vjd'^Vi = -dip + vA'v l + f i} (18) 

d[v, = 0. (19) 

This flow might be freely decaying, or it might for example 
be a steady (or statistically steady) flow sustained by a body 
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force. (Since we have omitted buoyancy forces, we cannot 
consider convection as such in this section, but we can study 
laminar or turbulent flows driven by a body force.) If we 
assume that d[fi = 0, since only the solenoidal part of the 
force drives the flow, then p satisfies the Poisson equation 

A'p = -(afojfl^i. (20) 

The Laplacian operator A' has a unique inverse A'^ 1 if, as 
we assume, the boundary conditions are periodic and the 
fields have zero mean. This inverse is also a self-adjoint op- 
erator. 

The linearized equations are 

(d' t + Vjd'j)5vi + (Svjd'j — av x d' y )vi + av x 5i2 

= -d[8p + aSud'yP + u(A'Svi - 2ad' x d' y Vi) , 

d-Svi - ad' y v x = 0, (21) 

where 8vi is the velocity perturbation induced at first order 
by the shear, and 8p is the accompanying pressure pertur- 
bation. In general the linearized equations must be solved 
numerically. 



3.2 Asymptotic analysis for high frequencies 

We now consider a limit in which the shear is oscillatory 
with high frequency. By this we mean that the timescale of 
the oscillations is small compared to the circulation period 
of the streamlines of the basic flow and the viscous timescale. 
For multiscale flows such as fully developed turbulence, the 
approximation considered here applies when the oscillation 
period is short compared to the convective timescale of the 
smallest eddies. Alternatively, it can be viewed as determin- 
ing the response of those eddies for which the convective 
timescale is long compared to the oscillation period. 

We use the method of multiple scales and introduce a 
fast time variable T' = t' /e, where e -C 1 is a small parame- 
ter that characterizes the ratio of timescales. The rapidity of 
the shear is expressed by rewriting a i-> a(T') and a e~ a, 
where the new meaning of a is da/dT' . The basic flow may 
vary with the ordinary time variable t'. We then expand 

5vi = Svio + eSvn + ■ ■ ■ , (22) 

Sp = e- 1 (Sp + eSp 1 + ■■■), (23) 

in asymptotic series, where the quantities on the right-hand 
side depend on (x' ,t' ,T'). The Sp expansion is indexed in 
this way for reasons that will become clear. Since the ba- 
sic flow does not depend on e, which is a property of the 
perturbations only, it is not expanded. 
At leading order we find 

d' T Sv i0 + av x Si2 = -d'iSp , (24) 

d'iSvio - ad' y v x = 0. (25) 

Roughly speaking, these equations describe an elastic re- 
sponse in which (from the second component of equation l24[) 
5v y o w —av x , leading to a shear stress —(v x 8v y o) « a(v x ) 
proportional to the shear and in phase with it. The real- 
ity is somewhat more complicated because of the pressure 



terms and the constraint of incompressibility. The pressure 
perturbation satisfies a Poisson equation, 

A'Spo = — 2ad' y v x , (26) 

obtained by eliminating Sv from the above two equations, 
and so 

d' T Svio = 2ad' i d y A'~ 1 v x — av x 8 i2 . (27) 
The linearized shear stress at this order is 

- 5R xy0 = ~(v x 5v y0 + v y 5v x o) (28) 

and satisfies 

d' T {-5R X yQ) = a(vl - 2(v x d'y + v y d' x )d' y A'~ 1 v x ). (29) 
In terms of the tensor 

A ljk i = (vid'jd' k A'~ vi), (30) 

we have 

d' T {-SR X y ) = a(A lnl - 2A1221 - 2A 2 i 2 i). (31) 

This tensor has the symmetry property Aijki = Aikji, which 
follows immediately from its definition, and the contraction 
properties Aijkj = Aijkk = 0, which follow from d[vi = 0. 
The further symmetry property Aijki = Aijki follows from 
integration by parts, if the averaging operation includes a 
spatial average over a periodic cell. The quantity multiplying 
a on the right-hand side of equation (|3ip can be regarded 
as the effective elastic modulus (shear modulus) of the flow 
(divided by the density of the fluid). 
At the next order we find 

d' T 5va + (d' t + Vjd'j)Svio + {5v j0 d'j - av x d' y )vi 

= -d'iSpi + aSud'yP + v(A'8v i0 — 2ad' x d'yVi), 

d'iSvn = 0. (32) 
We then obtain the Poisson equation 

A'Sp! = -a(d' t + Vjdj)d y v x - 2(d' i Vj)d' j 8v i0 
+ a(d[v x )d'yVi + ad x d' y p + avA'd' y v x , 

which can be inverted in principle. The linearized shear 
stress at this order is 

— 8R xy i — —(v x 8v yl + v y 8v x i) (33) 

and satisfies 

d' T (-SR xyl ) = {v x (d' t + Vjd'^Svyo + v y (d' t + Vjd'^Svxo 
+ v x (8v j0 d' j - av x d' y )vy + Vy(5vj d'j — av x d' y )v x 
+ (v x d' y + v y d x )5pi - avyd'yP 

— vv x {A'8v y o — 2ad' x d y Vy) 

- vv y (A'8v x0 - 2ad' x d' y v x )) . 

We substitute for p and Spi from the Poisson equations that 
they satisfy, using the fact that the inverse Laplacian is self- 
adjoint, and integrating by parts in various places: 

d' T {-5R xv i) = (v x (d' t + Vjd'j)8v y o + v y (d' t + Vjd'j)5v x o 
+ v x (5vjod'j — av x d'y)v y + v y (5vjod'j — av x d' y )v x 

- A'~ 1 (d' y v x - d' x v y )[-a(d' t + Vjdj)d'yV x 

— 2(d' i Vj)d' j Sv i o + a(d'iV x )d'yVi + ad' x d' y p + avA'd' y v x ] 

— a(d' i Vj)(d'jVi)dyA'~ l v y — uv x (A'Sv y0 - 2ad' x d' y v y ) 

- vv y {A'5v x0 - 2ad' x d' y v x )). 
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A further replacement of p is required: 

d' T (-5R xyl ) = (v x (d' t + Vjd'j)5v y0 + v y (d' t + Vjd'j)6v x0 
+ v x (5vjod'j — av x d'y)v y + v y (8vjodj — av x d' y )v x 

- A'^ 1 (8' y v x + 8' x v y )[-a(d' t + Vjd'j)d y v x 

— 2(d'iVj)d'jSvio + a(d'iV x )d' y Vi + avA'd' y v x ] 
+ aA'~ 2 d' x d y (d' y v x + d' x v y )(d' i v j )(8' :j v i ) 

— a(d' i Vj)(djV i )d y A'~ 1 v y — uv x (A'Sv y0 - 2ad' x d' y v y ) 

- vv y (A'5v x0 - 2ad' x d' y v x )) . 

Since we have an expression for d' T Svio, we consider 
&T{-5R xyl ) = {v x (d' t + Vjd'j)d' T 8v yQ + v y (8' t + Vjdj)d' T Sv x o 
+ v x [(d' T 8vj )d'j - av x dy]v v + v y [(d' T 8vjo)d'j — av x 8' y ]v x 

- A'~ x (tf y v x + d' x v y )[-a{d' t + Vjd'^d'yVx 

- 2(d'iVj)djdTdvio + a{d'iV x )d' y Vi + auA'd' y v x ] 
+ tiA'~ 2 d' x d' y (d y v x + ^%)(5-u f )(^«i) 

- a(d'iVj)(djVi)d' y A'~ v y ~ vv x (A' 'd' T 8v yQ — 2ad' x d' y v y ) 

- vv y (A' 'd' T 8v x0 ~ 2ad' x d' y v x )) . 
Substituting for d' T 5vio, we obtain 

dT{SR xyl ) = a(v x (d' t + Vjd' j )(2d' y d' y A'~ 1 v x - v x ) 
+ v y (d' t + v ] d' J )(2d' x d' y A'- 1 v x ) 
+ Vxl^d'jd'yA'" 1 v x )d'j — 2v x d' y ]v y 
+ Vyl^d'jd'yA'" 1 v^d'j - 2v x d' y \v x 

- A'^ 1 (d' y v x + d' x v v )[-(d' t +v J d'j)d' y v x 

- 2(d' i Vj)d' :j {2d'id l y A'" 1 v x - v x S i2 ) + (d'iV x )d' y Vi 

+ vA'd' y v x ] + A'~ 2 d' x d' y (d' y v x + d' x v y )(d' i Vj)(d' j Vi) 

- (d' i v J )(d' j Vi)d y A'~ 1 v y 

- vv x (2d' y d' y v x - A'v x - 2d' x d' y v y )). 

This can be rearranged as follows, after integration by parts: 
d'T(-5R xy i) = a((d' t v x )(-v x + d' y d' y A'~ 1 v x + d' x d' y A'~ 1 v y ) 

+ v x {vjd'j){2d'yd' y A'~ 1 v x - v x ) 

+ v y {vjd'j){2d' x d' y A'~ 1 v x ) 

+ Vx^d'jd'yA'^ 1 v x )d'j - 2v x d y ]v y 

+ Vyl^d'jd'yA'' 1 v x )d'j - 2v x d' y ]v x 

- A'^ 1 (d' y v x + d' x v y )[-(vjd'j)d' y v x 

- 2(d[v j )d' J {2d' l d' y A'-\ x ) + 3(8^)8'^] 
+ A'~ 2 d' x d' y (d' y v x + d^v^id'iV^idjVi) 

- (d' l v j )(d' j v i )d' y A'~ 1 v y 

- vv x (d' y d' y v x - A'v x - 3d' x d y v y )}. 
We define the tensors 

Bijki = {(d' t Vi)d' j d' k A'~ 1 vi), 



Cijki = -v(vid'jd' k vi), 



Dijki = (viVjd' k vi}, 



Dijkimn = (viVjd k did m A v n ) 



(34) 
(35) 
(36) 
(37) 



Eijki = (v m (d' m d' n A' 1 d' i v j )d' n A' 1 d' k vi), (39) 

which satisfy B t jki = Bikji, Bij k j = 0, Cijki = Cikji = Cijki, 
Cijkj = and various other identities for D and E. We then 
find (after integration by parts) 

8't (— SR xy i) = —a(B\jj\ — B1221 — -B1122 

— C1221 + Cijji + 3C1122 

— 2Dljj221 — 2_D2jjl21 — 3Dijj221 — 3Dijj212 

— Dijijl221 — Dijij\212 + Dijij22 
+ 4i?2121 + 4i?2112)- 



3.3 Interpretation 

The two results obtained so far can be written in the forms 
d'ASRxyo) = aGo, (40) 



d' T {—8R xy i) = —aQi, 



(41) 



where Go and Gi are two coefficients, each of which could be 
positive, negative or zero. Given that this a linear analysis, 
we may consider a complex shear a tx exp(— iwi) with angu- 
lar frequency cj = 0{e~ 1 ) and deduce that the shear stress 
in the limit of high frequency is 



5R X 



Go 



(42) 



The term Go represents an ideal elastic response, while the 
term Gi represents an imperfection of the elasticity associ- 
ated with dissipation. For comparison, the shear stress of a 
viscous fluid is va = — iuji/a. Therefore the effective kine- 
matic viscosity of the flow at high frequ e ncies is Gi/oj 2 . 

The calculation of iGoodman fc Ohl l| 19971 ) corresponds 
only to the leading order of the above expansion. They do 
not mention the elastic stress but focus on the dissipation 
rate at leading order in a periodic strain, given by their 
equation (25). As they point out, this expression evaluates 
to zero; this is consistent with our analysis because there is 
no dissipation associated with a perfect elastic stress. They 
obtain a non-zero result at this order by manipulating the 
apparent singularity in their expression at zero frequency. 
While this procedure may produce a meaningful result, it is 
not really justified because the preceding steps have assumed 
a separation of timescales between the tide and the convec- 
tion; the zero-frequency pole signals the breakdown of their 
approximation scheme and its resolution is not straightfor- 
ward. In contrast, our calculation of Gi and the effective 
viscosity in the high-frequency limit is based on a system- 
atic asymptotic expansion. 



3.4 Evaluation for statistically isotropic flows 

Although convective flows are naturally anisotropic, analyt- 
ical progress is easiest when the flow is assumed to have 
isotropic properties when spatially averaged. In this case 
Aijki must be an isotropic tensor, i.e. 



Aijki — AiSijSki + A25ikSji + A?,8u8 



(43) 



klmnpq = (viVjd'kd'td^d'nd'pA' v q ), (38) 



The symmetry and contraction properties then imply A2 = 
Ai and A3 = — (d + l)Ai, where d is the number of spatial 
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dimensions (of course d = 3, but the case d = 2 is also of at 
least theoretical interest), and so 

2K 

Aijkl = d(d -l){d + 2) K d + 1 ) g<tg J* ~ SijSkl ~ ^J'L ( 44 ) 
where the mean kinetic energy density if is given by 

2K = {viVi) = = -d(d -l)(d+ 2)Ai. (45) 
In this isotropic case we then obtain 



_ 2(d-2)(d+l) 
W0 d(d-l)(d + 2) ' 



(46) 



The effective elastic modulus is positive for a three- 
dimensional flow but vanishes in two dimensions under the 
assumption of isotropy. 

In the isotropic case it can be shown (see Appendix [Bjl 
that the triple-correlation tensor D vanishes identically. E 
also vanishes when the identity Eij^i = —Eunj (which fol- 
lows from integration by parts) is combined with the general 
form of a fourth-rank isotropic tensor. The tensors B and C 
have the form 



Bijti — B[(d+ l)Su8jk — SijSki — SikSji], 



(47) 



Cijki = C[(d + l)5 a 8 jk — SijS kt - S ik Sji]. (48) 
Furthermore, the energy equation of the basic flow, 

= v(viA'vi) + {fiVi}, (49) 

implies 



B, 



~\~ Cijji — (fiVi) — Py 



(50) 



the power input per unit volume (or area), and so 

d(d-l){d + 2)(B + C)=P (51) 
for isotropic statistics. In this case 



& = B(d 2 -2) + C(d 2 -6). 



(52) 



This corresponds to a viscous response, with the effective 
viscosity being inversely proportional to frequency-squared, 
for high-frequency oscillatory shear. The effective viscosity 
coefficient is Qi/ui 2 . In the case B = 0, when the flow is 
maintained steadily against viscous dissipation by the body 
force, the effective viscosity is positive in three dimensions 
and negative in two dimensions. In the case P — 0, when the 
flow is freely decaying (-B = C > 0), the effective viscosity 
is negative in either three or two dimensions. 



3.5 Evaluation for ABC flows 

A widely studied class of incompressible fluid flows in a 
periodic domain is provided by the AB C flow (named 
after Arnol'd, Beltram i and Childress; lArnol'dl 1 19651 ; 
iGallowav fe Frischll 19871 ) 



(A sin kz' + C cos ky ' ^ 
B sin kx' + A cos kz' 
C sin ky' + B cos kx' i 



(53) 



in a cube of length 2-k jk. This velocity field has the Beltrami 
property V' X v = kv, so nonlinearity is absent; v ■ V't) is 
balanced by a pressure gradient. If the flow is unforced, A, 
B and C decay proportionally to exp( — vk 2 t). Alternatively 



the flow can be maintained against dissipation by supplying 
a body force / = vk 2 v. The most widely studied example 
has A — B — C. 

The response coefficients are easily evaluated as 



G Q \(A 2 -C 2 ), 
Gi = \A{A + vk 2 A). 



(54) 



(55) 



Therefore the elasticity can be positive, negative or zero de- 
pending on the anisotropy of the flow and its orientation 
relative to the shear. The effective viscosity Q\/ui 2 at lead- 
ing order vanishes for a freely decaying flow but is posi- 
tive (and inversely proportional to frequency-squared) for a 
forced flow, assuming that A 0. 

In fact, probably all analytical examples of three- 
dimensional fluid flows lack genuine nonlinearity, in the sense 
that, if they are expanded in a Fourier basis with wavenum- 
bers fe, there are no non-empty triads of interacting com- 
ponents satisfying fci + k-2 + ks = 0. Such triads tend to 
produce a cascade of energy to larger wavenumbers in the 
manner of hydrodynamic turbulence. If there are no non- 
empty triads then the triple-correlation tensors D and E 
again vanish, and the only contributions to Q\ come from 
the time-dependence of the flow (the tensor B) and the vis- 
cous terms (the tensor C). Both of these effects may be small 
if the viscosity is small. 

The ABC flow is stable only for sufficiently small 
Reynolds number. More realistically, in a typical flow at high 
Reynolds number there will be a turbulent cascade involving 
strong triad interactions, meaning that the tensors D and 
E may be significant (although apparently not in isotropic 
turbulence). The tensors B and C will also be enhanced by 
the turbulent cascade. Numerical simulations are required 
to access this regime. 

We note that the 'eddy viscosity' of the A = B — C 
flo w has been calculated, as a fu nction of Reynolds number, 
bv lWirth. Gama fe Frischl (|l995l ); see also references therein 
for related studies. However, the nature of their calculation 
is different; using multiscale techniques, they determine the 
behaviour of very slow large-scale deformations of the cel- 
lular flow, and deduce that a large-scale instability occurs 
through the appearance of a negative eddy viscosity. In con- 
trast, our analysis determines the response of a flow to an 
imposed large-scale deformation of high frequency. 



4 DIRECT NUMERICAL SIMULATIONS OF 
CONVECTION IN AN OSCILLATORY 
SHEARING BOX 

4.1 Numerical setup 

We have implemented the oscillatory shearing box in the 
SNOOPY code, a 3D spectral code solving the equations of 
incompressible or Boussines q (magneto)hydrodynam i cs us- 
ing a Fourier represent ation (|Lesur fe Longa rctti 2005. [20071 ; 
iLesur fe Ogiiviel |2010| V We assume that an external force 
creates an oscillatory shear with d = — 5'cos(Slt), where S is 
the maximum shear rate and Q is the tidal frequency. As in 
previous sections, the shearing motion is in the y direction, 
with a linear dependence on x, i.e. — S cos(£lt)x e y , and the 
force that drives it is fiS 1 sm(Qt)x e y . 
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In addition to this imposed shearing motion, convection 
is driven by applying uniform gravity and an unstable en- 
tropy gradient, within the Boussinesq approximation. (No 
additional body force is required to drive the motion.) By 
setting N 2 = — 1, we adopt a unit of time related to the 
unstable stratification. We test two different configurations: 
convection in the shearwise direction (e = e x ) and con- 
vection in the spanwise direction (e = e z , where e is the 
direction of stratification defined in Section [3J. 

The aspect ratio of the periodic box is adjusted accord- 
ing to the direction of stratification, with L x x L y x L z = 
2x1x1 when e = e x and 1x1x2 when e — e z . Hav- 
ing a box elongated in the direction of stratification allows 
'elevator modes' to break up more easily through secondary 
instabilities; otherwise, these tend t o dominate the convec - 
tion at moderate Rayleigh numbers l|Lesur fe Qgilviell2010l ). 
The resolution is 128 collocation points per unit of length. 

To control dissipation on small scales, we introduce an 
explicit kinematic viscosity v and thermal diffusivity n, with 
Prandtl number v/k = 1 for simplicity. The value of the 
diffusion coefficients is set according to the Rayleigh number 
Ra = \N 2 \L y /vn. In our setup, convection starts when Ra > 
(2tt) 4 « 1559. In the following we will consider simulations 
exhibiting fully turbulent convection, with Ra = 4 x 10 6 . 

In order to avoid any artefact of the initial conditions, 
we initiate our simulations with noise at the largest scales 
and we let turbulence evolve without any shear for 100 
turnover times (i.e. from t = — 100|7V| _1 to t = 0). The 
spectrum of the turbulence we obtain and a typical snap- 
shot are shown in Fig. [T] Once a quasi-stationary turbulent 
state is reached, we switch on a weak oscillatory shear and 
start measuring the Reynolds stress —R xy (t) = —{v x v y ), 
where (•) denotes a volume average over the box. Such a 
simulation has to be continued for an integration time AT 
of hundreds of turnover times in order to reduce the impact 
of turbulent noise on the measurements. 

4.2 Turbulent viscosity: definitions and simple 
models 

Traditionally, the turbulent viscosity is associated with a 
simple closure formula relating the Reynolds stress to the 
rate of strain; in our system, 

— R xy = v t a. (56) 

This expression is based on the assumptions that the re- 
lationship between stress and rate of strain is linear and 
instantaneous. In the case of rapid oscillatory shear, the as- 
sumption of linearity is reasonable provided that the shear 
is of sufficiently small amplitude. However, the assumption 
of instantaneity is generally not valid and we should allow 
the Reynolds stress to be linearly related to a(t) through an 
integral over its previous history. In the Fourier domain, this 
relationship (a convolution) reduces to a multiplication: 

- R xy {uj) = z/t(w)d(w), (57) 

where ~ denotes the Fourier transform in time. In this ex- 
pression Vt(u) is a complex function of frequency, rather 
than a real constant, indicating that the Reynolds stress 
and the rate of strain are generally out of phase, and that 
the relationship is frequency-dependent. Since the relation- 
ship involves real-valued functions in the temporal domain, 




k/2n 



Figure 1. Simulation snapshot of thermal fluctuations (top) and 
energy spectra (bottom) at Ra = 4 X 10 6 in the case without 
any shear. We show the kinetic energy spectrum (E^) and the 
temperature fluctuation spectrum (6fe). 

vt has the Hermitian symmetry ft(— u>) = [i/t(w)]* for real 
uj. With our choice of a = — Scos(fit), we have 

R xy (u}) = v t {uj)SiT[6{w - fi) + 5{u + fi)]. (58) 

This makes it possible to measure v x (fi) from the time series 
R xy (t). In reality, the delta function is replaced by a peak of 
finite height and non-zero width because of the finite inte- 
gration time; alternatively, it becomes a Kronecker delta in 
a discrete Fourier transform. Furthermore, noise is present 
because there are turbulent fluctuations in R xy even in the 
absence of shear. The integration time must be sufficiently 
long to allow the signal to be detected above the noise. 

A simple closure model for the Reynolds stress that can 
be compared with the results of numerical simulations is a 
viscoelastic model, which reflects the fact that the turbulent 
stress cannot instantaneously produce a viscous response to 
a time-dependent shear rate. [Viscoelastic models for magne- 
tohydrodyn amic turbulenc e in astrophysical discs have been 
discussed bv lQgilviell|200ll ).] We suppose that the convective 
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turbulence contains a number of viscoelastic components, 
having effective elastic moduli Cj and relaxation times 1 j^j , 
giving rise to the complex turbulent viscosity 

d 



ft(w) = ^2 



11 



(59) 



In this expression, Cj and jj are real, but we do not require 
Cj > 0, allowing for the possibility that negative Re[ft(w)] 
may appear for some ranges of frequency. For example, the 
case of constant shear rate (uj = 0) is known to produce 
negative viscosity for low values of Ra in protoplanetary 
disc convection, which involves rotation as well as shear 
l|Lesur fc QgilvitfeOld ). However, in order to have a physical 
relaxation for all possible input signals, we require jj > 0, 
indicating that the system we describe is stable for any time- 
dependent shear. In other words, since Vt(ui) in equation (|57[) 
derives from a causal integral relationship, it should be an- 
alytic in the upper half-plane. 

This model shares several properties with the analytical 
results derived in Section[3] In the high-frequency limit cj — > 
oo, the model gives 



(60) 



This behaviour is equivalent to equation (|42p obtained 
for the high-frequency response of an arbitrary flow, if 
Go — Y2j c i (t ne high-frequency elastic modulus) and Q\ — 
X]j c i7i- (O n the other hand, the viscosity obtained in the 
low- frequency limit w — ► is J^. Cj/^j.) As we will see later, 
the same type of dependency is found in numerical simula- 
tions. 

In this paper, we will assume that only two viscoelastic 
components are present. This can be seen as a computational 
limitation, since numerical simulations cannot probe very 
high frequencies which are usually associated with scales be- 
low the grid size. In the following, we will therefore consider 
the simplified model 

ci c 2 



+ 

71 — iuj 72 — iui 



(61) 



4.3 Measuring the turbulent viscosity 



A typical example of the time series of the Reynolds stress 
and its Fourier transform are shown in Fig. [2] In order to re- 
duce the aliasing of low frequencies into the high-frequency 
domain due to the finite integration time, we have applied 
a Hanning window to the time series before computing 
the Fourier transform. As is evident from the time series, 
the turbulent convection produces large fluctuations in the 
Reynolds stress which dominate the response to the oscilla- 
tory shear. Nevertheless, it is possible to extract useful infor- 
mation by looking at the temporal spectrum of the stress. In 
particular, a spike is clearly visible at ui — Q, which indicates 
that the turbulent flow is producing a detectable response 
to our forcing. With the complete Fourier transform of the 
Reynolds stress, it is possible to measure the turbulent vis- 
cosity at the forcing frequency using equation (|58p . More- 
over, evaluating the noise level in the vicinity of the spike 
allows us to estimate the uncertainty in the measurement of 
the turbulent viscosity. 

It should be noted that the signal-to-noise ratio (de- 
fined as the ratio of the amplitude of the spectral spike to 



the amplitude of the surrounding noise; see Fig. [2] right 
panel) should decay as AT~ 1//2 . Therefore, if the oscillatory 
signal (the spectral spike) is too weak to be detected in the 
Fourier transform, one can in principle increase the integra- 
tion time to reduce the signal-to-noise ratio, but this is a 
computationally expensive procedure. 



4.4 Numerical results 

The results of the simulations performed in this work are 
shown in Table [T] They are also plotted in Fig. [3] for shear- 
wise convection and in Fig. [4] for spanwise convection. De- 
spite the present of significant noise in some of the results, 
several deductions can be made from these simulations: 

• At sufficiently high frequency, we find Re(«/ t ) < in 
every case. Note, however, that noise dominates many high- 
frequency measurements of this quantity. To substantiate 
the high-frequency behaviour, we have performed a simula- 
tion at lower Rayleigh number (Ra = 10 5 ) for more than 
3000 turnover times in order to reduce the effect of turbu- 
lent noise. This simulation exhibits clearly Re(f t ) < at 
Q = 64, indicating that the trend observed in the other con- 
figurations is real. 

• Im(ft) > 0. This conclusion is rather strong since the 
noise level is smaller than the measured value of Im(f t ) in 
most cases. It means that the effective elasticity is positive. 

• Larger values of \i/ t \ are found when the stratification 
is in the x direction (shearwise convection). This suggests 
that the full turbulent viscosity tensor is anisotropic. 

• Changing 5* does not change significantly the measured 
turbulent viscosity, suggesting that the numerical simula- 
tions are in the linear regime assumed in writing down equa- 
tion (|57|l . However, the results obtained with the smaller 
value of S are subject to significant uncertainty. 

• The asymptotic behaviour in the high-frequency limit 
is Re(^t) oc cj -2 and Im(f t ) oc cj" 1 . 



We have also fitted the simple closure model (|61[) to 
the numerical results for shearwise and spanwise convection. 
The best fits are shown as green curves in Figs [3] and [4] and 
the coefficients we obtained are shown in Table [2] The ma- 
jor viscoelastic component has an effective elastic modulus 
ci that is comparable to, but somewhat less than, the mean 
kinetic energy density K (ci. eq. 1461 which suggests a ratio 
of 4/15 in the isotropic case) and a relaxation rate compara- 
ble to the nominal convective turnover rate \N\, which is 1 
in our system of units. Note that the contribution to Re(ft) 
from the first component is (ci/7i)[l + (w/71) 2 ] -1 , which 
is compatible with the form of the frequency-dependent vis- 
cosity used bv lQgilvie fcTii] \200H ). The suggested presence 
of a second viscoelastic component with a negative C2 and a 
faster relaxation indicates that, at high frequency, the real 
part of the viscosity changes sign, as is observed in the nu- 
merical results. With more data we might be led to introduce 
a broader spectrum of viscoelastic components. 

In order to identify which spatial scales contribute to 
the turbulent viscosity, we have computed the spatial 'spec- 
trum' of the turbulent viscosity. This is derived from the 
instantaneous spatial spectrum of the Reynolds stress, 



R xy (k,t) = 2Re[# xc '0*j, 



(62) 
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Figure 2. Time history of the Reynolds stress (left) and Fourier transform of the Reynolds stress (right) with an oscillatory shear of 
frequency Q = 16 and amplitude S = 0.1. The periodic signal cannot be seen in the time series. The Fourier transform exhibits a peak 
at the excitation frequency u> = Q. 
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Table 1. Parameters and results of the simulations discussed in this paper. The error in ft is estimated based on the rms turbulent noise 
in the Fourier transform at the oscillation frequency. K is the mean kinetic energy density of the fluctuations. 



coefficient shearwise convection spanwise convection 

ci 0.0978 0.06 

71 1.07 1.62 
c 2 -0.0286 -0.0178 

72 6.15 24.25 

Table 2. Best-fit coefficients for the closure model H61I I. 

where i>i is the 3D instantaneous spatial Fourier transform 
of Vi and the overbar denotes a shell-integration proce- 
dure in spectral space, such that R xy (t) — R xy (k,t) dk. 



R xy (k,t) is then used in a relation similar to ([58} which 
defines the turbulent viscosity spectrum ut(k,u>): 

R xy {k,uj) = v t {k,uj)S 7r[<5(o; - O) + 5{u + O)] . (63) 

We have applied this procedure to the simulation with 
Q = 16, S = 0.1 and with stratification in the z direction. 
The corresponding spectra are shown in Fig. [S] As can be 
seen, much of the turbulent viscosity (both real and imagi- 
nary) is due to the largest scales of the system. However, we 
also observe a component with Re(^ t ) > at 'small' scales 
(k/2n ~ 10), indicating that there may be an important 
contribution from scales whose turnover time is of the or- 
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Figure 3. Turbulent viscosity of shearwise convection, versus the angular oscillation frequency in units of \N\. The numerical mea- 
surements are shown in blue (open circles) and the best fit of the closure model d 6 1 1) is shown in green (triangles). Negative values are 
connected by a dashed line. Uncertainties due to turbulent noise are shown as a shaded region. 
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Figure 4. Turbulent viscosity of spanwise convection, versus the angular oscillation frequency in units of |7V|. The numerical measure- 
ments are shown in blue (open circles; S = 0.1) and red (crosses; S = 0.05) and the best fit of the closure model l|61|l is shown in green 
(triangles). Negative values are connected by a dashed line. Uncertainties due to turbulent noise are shown as a shaded region. 



der of the tidal frequency. It should be stressed that these 
spectra are strongly polluted by turbulent noise, especially 
at low wavenumbers. Therefore, this result should be seen 
as a plausible trend. Longer simulations having a larger res- 
olution should be used to confirm this finding. 



5 SUMMARY AND DISCUSSION 

In this paper we have studied the interaction between tides 
and convection in astrophysical bodies by analysing the ef- 
fect of a homogeneous oscillatory shear on a fluid flow. This 
model can be taken to represent the interaction between 
a large-scale periodic tidal deformation and a smaller-scale 
convective motion. We first considered analytically the limit 
in which the shear is of low amplitude and the oscillation pe- 
riod is short compared to the timescales of the unperturbed 
flow. In this limit there is a viscoelastic response and we ob- 
tained expressions for the effective elastic modulus and vis- 
cosity coefficient. The effective viscosity is inversely propor- 
tional to the square of the oscillation frequency, with a co- 



efficient that can be positive, negative or zero depending on 
the properties of the unperturbed flow. We also carried out 
direct numerical simulations of Boussinesq convection in an 
oscillatory shearing box and measured the time-dependent 
Reynolds stress. The results indicate that the effective vis- 
cosity falls rapidly as the oscillation frequency is increased, 
attaining small negative values in the cases we have exam- 
ined. 

Our methods and findings diffe r signi f icantl y from those 
of other authors. The hypothesis of lZahnl l|l966l ) that the ef- 
fective viscosity of large eddies is reduced by only a linear 
factor at high frequencies is incompatible with our analyti- 
cal and numerical results. We find a quadratic reduction fac- 
tor at high frequencies, similarly tolGoldreich fc Nicholso3 
(|l977h and iGoldreich fc Keelevf (|l977t ), but not necessar- 
ily for the sam e reason. Our analytica l study, which goes 
beyond that of iGoodman fc Ohl (| 19971 ), implies that large 
eddies generically provide a quadrat ically reduced viscosity 
[which lGoldreich fc Nicholson! (11977T ) considered as an upper 
bound], but with a coefficient that can be positive, negative 
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Figure 5. Spatial spectrum of the turbulent viscosity for the simulation of spanwise convection with f2 = 16 and S = 0.1. The real part 
has two components: a negative contribution from large scales and a positive contribution from smaller scale (k/2n ~ 10). 



or zero. Our numerical results also indicate a quadratically 
reduced viscosity, with a tendency towards negative values 
at high frequencies in the cases we have examined, an effect 
due to the largest scales in the turbulent flow. 

The appearance of a negative effective viscosity may be 
alarming. However, there is no reason in principle why the 
effective viscosity of a convective fluid need be positive. The 
consequence of a negative effective viscosity would be that 
tidal evolution proceeds in the reverse direction to that usu- 
ally assumed. Angular momentum would be transferred in 
an 'anti-frictional' direction, from the less rapidly rotating 
component to the more rapidly rotating one. For example, if 
a small satellite orbits a body with a negative effective vis- 
cosity, and if the larger body spins more slowly than the or- 
bit, then the orbit would expand and the larger body would 
spin down. Energy would be added to the spin-orbit system. 
(Similarly, orbital eccentricity could be caused to increase in 
situations where circularization would usually be expected.) 
While this effect might be called tidal 'anti-dissipation', it 
does not contradict the laws of thermodynamics. Within lim- 
its, work can indeed be done by the convection on the tidal 
flow; the energy comes from the buoyancy forces that drive 
the convection, and ultimately from nuclear or gravitational 
energy. Since we are considering a regime in which the tidal 
strain is small and the tidal frequency is high (implying that 
the negative effective viscosity is much smaller in magnitude 
than the values estimated from mixing-length theory), only 
a very small fraction of the energy budget of the star would 
be diverted into the spin-orbit system. The consequences 
could still be important, because the nuclear energy content 
of the Sun, for example, is about five orders of magnitude 
larger than the orbital energy of (say) a solar-type binary 
star with an orbital period of ten days. 

A well known example of a related phenomenon is the 
negative effective viscosity of convect ion at low or moderate 
Rayleigh number in an accretion disc |Lesur fc Ogilviell2010l . 
and references therein), which would lead to anti-frictional 
angular momentum transport and anti-diffusion of the sur- 
face density of the disc. This example also serves as a re- 
minder that the Coriolis force present in a rotating system 
can affect the transport of (angular) momentum; it would 
therefore be of interest to include rotation in the calcula- 



tions presented in this paper. Anti-frictional processes are 
also well known in the context of the Earth's atmosphere 
(e.g. lMcIntvre]|2000t ). 

There are several reasons why this reversed tidal evolu- 
tion may not, in fact, be found in nature. First, it may be 
that when a more realistic tidal velocity field is considered 
and the actual anisotropies of the convection and the Cori- 
olis force are taken into account, the net effective viscosity 
turns out to be positive. Second, in a solar-type star where 
the convective timescale becomes short near the surface, the 
net effective viscosity may be dominated by these more rapid 
convective elements and work out to be positive. Third, it 
may be that tidal dissipation is dominated by internal waves 
and that convection makes a relatively small contribution. 
Fourth, the negative effective viscosities that we find might 
be due to explicit (molecular) viscous effects which are to- 
tally negligible in astrophysical objects. Indeed, we have not 
found a negative turbulent viscosity that exceeds the ex- 
plicit (molecular) viscosity in magnitude. Nevertheless, it is 
of interest to note the theoretical possibility of reversed tidal 
evolution and to be alert to any observational evidence of 
such an effect. 

Penev and collaborators have attempted to address the 
issue of the reduction of the effective viscosity at high fre- 
quencies through simulations of convection in a deep layer. 
While their results suggest only a linear reduction factor 
(and a positive numerical value), there are many differences 
between their work and ours. Our numerical study involves 
a local model of convection, which is in some ways more 
limited but is also more controlled because a uniform un- 
stable stratification is imposed. We measure the effective 
viscosity directly rather than relying on a perturbative cal- 
culation or an indirect measurement. IPenev et alj (|2009h do 
not find good agreement between different methods, and 
their results depend on the shear amplitude. Our work cov- 
ers a wider range of oscillation frequencies and quantifies 
the uncertainty due to turbulent noise. It should also be 
pointed out that our work uses a very small forcing com- 
pared to Penev and collaborators. In particular, the rms ve- 
locities with and without fo rcing are identical in all our runs, 
whereas IPenev et all {2009) have forced velocities of the or- 
der of the rms velocity in their case of weak forcing, and 
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even 5 times larger in the case of strong forcing. We believe 
that such strong forcing could significantly affect the mea- 
surement of turbulent viscosity. Several tests that we have 
performed indicate that doubling the forcing amplitude (to 
S = 0.2) is already enough to modify our results. 

In our opinion more work is required to extend both 
analytical and numerical calculations of this type before the 
frequency-dependent effective viscosity of a realistic stellar 
or planetary convective zone can be reliably estimated, and 
the efficiency of tidal dissipation in astrophysical bodies can 
be better understood. 
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APPENDIX A: COMPOSITION OF AN ARBITRARY TRACELESS VELOCITY GRADIENT TENSOR 
FROM SIMPLE SHEARS 

We consider an arbitrary incompressible velocity gradient tensor, being a 3 x 3 matrix with vanishing trace. Each of the 
off-diagonal elements represents a simple shear, and therefore the off-diagonal part of the matrix is a linear combination of 
simple shears. To compose the diagonal elements, we first consider the matrix 

(Al) 

which is a linear combination of simple shears. By a rotation through n/4 about the z axis, we bring this matrix into the 
diagonal form 

/I 0\ 

-1 . (A2) 
\0 0/ 

In a similar way the matrix 

(\ \ 

(A3) 

can be constructed from simple shears. Using a linear combination of these two diagonal matrices, the three diagonal elements 
of the arbitrary velocity gradient tensor can be composed, subject to the constraint that the sum of the diagonal elements 
vanishes. 




APPENDIX B: VANISHING OF CERTAIN TRIPLE VELOCITY CORRELATIONS FOR 
STATISTICALLY ISOTROPIC FLOWS 

Let v(x) be a solenoidal vector field in d > 2 dimensions satisfying periodic boundary conditions, and let {■} be a spatial 
average over a periodic cell. (We omit primes on x and V in this Appendix.) We assume that v(x) has zero mean and 
is statistically isotropic so that any averages of tensor products of v and its derivatives are isotropic tensorsQ The inverse 
Laplacian operator A" 1 is well defined and self-adjoint for periodic fields of zero mean. 
We first consider the tensor 

D abcd = {v a v b d c v d ). (Bl) 
As an isotropic tensor, this must be a linear combination of products of Kronecker deltas: 

D a b c d = t\8ab&cd + t2&ac&bd + t3&ad&bc- (B2) 

The symmetry and contraction conditions D abc d = D bac d and D abcc = imply £2 = £3 and dti + £ 2 + £3 = 0. Thus 

Dabcd = t4.\d(8 a c8bd + S a dSbc) - 25 ab 8 cd ]. (B3) 

However, we also have 

D a bba = (VaV b d b Va) = (d b ( \v a v a v h ) ) = = £ 4 d(d - l)(d + 2), (B4) 

and so 

Dabcd = 0. (B5) 

We next consider 

Dabcdef = {v a V b d c dddeA~ 1 Vf). (B6) 

Isotropy and symmetry imply 

Dabcdef = £5 {8 a b8cd8ef + 8ab&ce8df + SabScfSde) 

+ t6(5 ac 8 b d8 e f + SacSbeSdf + SadSbcSef + S a dSb e S c f + SaeSbc^df + 5ac&bd&cf) 
+ tr(5acSbfSde + 8 a d5bf5ce + S a eS b f5 c d + S a f5 b c5de + S a fSbdS C e + SafSbeScd)- 

(There are three types of term here: those in which the two undifferentiated vs are paired, those in which each such v is paired 



1 It is debatable whether this assumption is compatible with the periodic boundary conditions. While the cubic symmetry imposed by 
the boundary conditions of a periodic cube are compatible with isotropy for tensors of second rank, this is not true for higher ranks. One 
could argue that the flow can be nearly statistically isotropic if it is dominated by scales smaller than the box. 
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with a d, and those in which one of them is paired with the differentiated v. Symmetry demands that the coefficients of terms 
of the same type are equal.) 

We require the contraction D abcdee (and those related by symmetry) to vanish. Thus 

(d + 2)ts5ab6 C d + te(d + 2)(5 ac Sbd + 8 ad 8 bc ) + 2t 7 (8 a bScd + S ac S bd + 8 a d8bc) = 0, (B7) 
which implies t& = t% and (d + 2)t 5 = — 2tv, and so 

Dabcdef = ts[(d + 2)(SacSbfS de + 5 a dSbfS C e + 8 a eSbf5 c d + 5 a fSb c Sde + 8 a fSbd5ce + 8 a fSb e 5 c d) 

— 2(S a b5cd5 e f + SabSceSdf + 5 a bS c fSde + 5 ac 8bd3 e f + 5 ac 8b e 5df + SadSbcSef + S ad 5beS c f + 5 ae 8 b cSdf + S ae 8bd5 c f)}. 

The contraction D a bcddf (and those related by symmetry) produces a Laplacian, so D abc ddf = D a b c f, which we have 
already shown to vanish. Thus 

(d + 4)t 8 [d(5 ac 6 b f + S a fSbc) - 25 ab S cf ] = 0, (B8) 

which implies tg = and so 

Dabcdef = 0. (B9) 

The last tensor to consider is 

D a bcdef gh = {v a v b dcd d ded f dg A~ 2 « h ) . (BIO) 
We give an abbreviated argument. Isotropy and symmetry imply 

Dabcdef gh = t$(5 a b&cd&ef & gh + 3 a bScdSegSfh + ■■■)+ tlo{Sac5bd5efSgh + 5 a c5bd8eh8fg + • • • ) 
+ tll(S a hSbc8dfSeg + &ah&bc&dg&ef + •'')• 

Requiring the contraction Dabcdef gg to vanish implies tg = tio and (d + 4)tg — —tn. The contraction D a bccef g h should produce 
D abefgh, which we have already shown to vanish. This is of the form given above, with t$ — (d + 4)tg + 2tio, te — (d + 6)£io 
and tj = 2tio + (d + 4)tn. Therefore tg — tio = tn = and so D a bcdef g h = 0. 
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